%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Half-Wave Semiconductor Nano-Laser %%%%
%%% v1 July 2023, S.H.Yun %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;

N = 20e5; %number of time steps
pumpwidth = 3; %ns
lambda_p = 1064; %[ns] pump wavelength
time = linspace(-pumpwidth*1, pumpwidth*1.2,N); %ns
pump = 2/sqrt(pi)* exp(-time.^2 * 4/ pumpwidth^2); %Gaussian pump profile
delt = (time(N)-time(1))/N;

cy = 0.92; %InGaAs semiconductor composition
Eg = 1.35 - 0.738*cy + 0.138*cy.^2; %Band gap energy in eV
m_e = 0.08-0.05*cy+0.017*cy^2; %effective mass of electron
m_h = 0.6-0.18*cy; %heavy hole mass, 0.12-0.1*y+0.03*y^2 for light hole
m = m_e*m_h /(m_e+m_h); %reduced carrier mass
h = 6.62e-34; %Planck's constant

M = 120; %Number of energy states in the model
energy = linspace(Eg, 1240/lambda_p, M); %spanning from the band edge to the pump wavelength

device = 'Device 3'; 
switch device
    case 'Device 1'
        L = 183; %Effective size of semiconductor particle [nm]
        thickness = 130; %Semiconductor thickness [nm]
        center = round(M*78/120); %position of the mode for 1200 nm
        Qfactor = 10;  %cold cavity Q factor
        Fp = 18; %Purcell factor at resonance peak
    case 'Device 2'
        L = 225; 
        thickness = 130; %Semiconductor thickness [nm]
        center = round(M*42/120); %position of the mode for 1350 nm
        Qfactor = 10;  
        Fp = 19; 
    case 'Device 3'
        L = 380; 
        thickness = 130; %Semiconductor thickness [nm]
        center = round(M*54/120); %position of the mode for 1305 nm
        Qfactor = 30;
        Fp = 8;
end

taus0 = 0.2; %ns, intrinsic fluoerscence decay time for free-space levels
Fp0 = Fp*taus0; %relative Purcell factor peak with reference to 1 ns
tauc0 = Qfactor*1240e-9/(2*pi*3e8)*1e9; %ns, cavity loss rate
taunr0 = 0.4e-4; %ns, non-radiative decay rate to lower levels
if delt/tauc0 > 0.2 || delt/taunr0 > 0.2, disp('Recommendation: Increase N'); end

V = L^2*thickness *1e-27; %semiconductor volume [m]
n0 = 1.5/M * 16*sqrt(2)/3 *pi /h^3 *(m_e*9.1e-31)^1.5* ...
    ((energy(end)-energy(1))*1.62e-19)^1.5 * V; %total number of states. m_e may be replaced by m
n0 = n0 * 0.5; %divided by two for only one spin (polarization) states
n_sat = n0 * ((0:M-1)'/M).^0.5; %The density of states of individual levels
pump0 = n0*M/1.5; %equal to the total number of carrier states

width = energy(center)/Qfactor /(energy(end)-energy(1)) * M; %the width of the single mode over the wavelength rage
shape = (width/2)^2 ./ (((1:M)' - center).^2 + (width/2)^2); %Lorentzian
x = (1:M)' - center;
B = width; %the coefficient in the stimulated emission term 
n = zeros(1,N); q=zeros(M,N); qtot=zeros(1,N);

K=3; %number of pump levels
px=linspace(0.8,8,K).^2; %Pump fluence range
pumppower= pump0 * px;
Pout=zeros(1,K); Ptot=zeros(1,K); fwhm=zeros(1,K); pk=zeros(1,K);
spectra=zeros(K,M); population=zeros(K,N);

tauc = tauc0*shape; %cavity lifetime parameter
tauc(delt./tauc > 1) = delt; %to avoid excessive change during one time step
Fp = Fp0*shape; %Frequency dependent Purcell factor profile
Fp(Fp < 1) = 1;
taus = taus0 ./ Fp; %spontaneous decay time
taus(delt./taus > 0.5) = 2*delt; %to ensure incremental changes per time step
taunr = taunr0 * ones(M,1);

for kk=1:K          %[ump fluence scan loop
    q=zeros(M,N);
    qtot=zeros(1,N);
    n =zeros(M,N);
    p = pumppower(kk) * pump;

    for ii=1:N-1
        taunr(M-1) = taunr0 ./ (n_sat(M-2) - n(M-2,ii) + 0.001) *n0;
        if n(M,ii)/taunr(M)*delt < n_sat(M-1) - n(M-1,ii)
            n(M,ii+1) = n(M,ii) + p(ii)*delt - n(M,ii)/taunr(M)*delt;
            n(M-1,ii+1) = n(M-1,ii) - B*n(M-1,ii)*q(M-1,ii)/taus(M-1)*delt - n(M-1,ii)/taus(M-1)*delt ...
                + n(M,ii)/taunr(M)*delt - n(M-1,ii)/taunr(M-1)*delt;
            q(M-1,ii+1) = q(M-1,ii) + B*n(M-1,ii)*q(M-1,ii)/taus(M-1)*delt + n(M-1,ii)/taus(M-1)*delt ...
                - q(M-1,ii)/tauc(M-1)*delt;
        else
            n(M,ii+1) = n(M,ii) + p(ii)*delt - (n_sat(M-1) - n(M-1,ii));
            n(M-1,ii+1) = n(M-1,ii) - B*n(M-1,ii)*q(M-1,ii)/taus(M-1)*delt - n(M-1,ii)/taus(M-1)*delt ...
                + (n_sat(M-1) - n(M-1,ii)) - n(M-1,ii)/taunr(M-1)*delt;
            q(M-1,ii+1) = q(M-1,ii) + B*n(M-1,ii)*q(M-1,ii)/taus(M-1)*delt + n(M-1,ii)/taus(M-1)*delt ...
                - q(M-1,ii)/tauc(M-1)*delt;
        end

        for mm = 2:M-2
            jj = M-mm;
            taunr(jj) = taunr0 ./ (n_sat(jj-1) - n(jj-1,ii) - n_sat(jj-1)*0.001) *n0;
            n(jj,ii+1) = n(jj,ii) - B*n(jj,ii)*q(jj,ii)/taus(jj)*delt - n(jj,ii)/taus(jj)*delt ...
                + n(jj+1,ii)/taunr(jj+1)*delt - n(jj,ii)/taunr(jj)*delt;
            q(jj,ii+1) = q(jj,ii) + B*n(jj,ii)*q(jj,ii)/taus(jj)*delt + n(jj,ii)/taus(jj)*delt ...
                - q(jj,ii)/tauc(jj)*delt;
        end

        qtot(ii) = sum(q(:,ii)); %total output photon over spectrum
    end

    spectra(kk,:) = sum(q,2)' ./ tauc'*delt; %output photon integrated over time
    population(kk,:) = sum(n,1); %total n integrated over spectrum

    figure(1); plot(energy, spectra(kk,:)); hold on;
    %figure(2); plot(time, qtot); hold on;
    figure(4); plot(energy, B*q(:,find(time>0,1,'first'))); hold on; %stimulated vs spontaneous ratio
    Qout(kk) = sum(qtot); %time integrated total intracavity photons
    Ptot(kk) = sum(spectra(kk,:)); %total extracavity photons
    Pout(kk) = max(spectra(kk,:)); %peak output power

    [mx, pk(kk)] = max(spectra(kk,:));
    y = spectra(kk,:)'/mx;
    ft = fittype('(w/2).^2 / ((x-c)^2 + (w/2)^2)'); %Lorentzian fitting function
    ff = fit(x,y, ft, 'StartPoint', [pk(kk)-center, 10]);
    z = (ff.w./2)^2 ./ ((x-ff.c).^2 + (ff.w/2)^2);
    fwhm(kk) = ff.w;
end

figure(1);
plot(energy, shape*max(spectra(K,:)) ,'k--');
xlim([0.77 1.165]);
xlabel('Energy (eV)'); ylabel('Output spectrum');
grid on; hold off;

%figure(2); xlabel('time'); ylabel('time-lapse intracavity power');
%grid on; hold off;

figure(3);
dE = (energy(end)-energy(1))/M*1e3; %meV
plot(energy, n_sat/dE, '--'); hold on;
plot(energy, n_sat/dE-n(:,find(time>0,1,'first'))/dE); hold on;
plot(energy, n_sat/dE-n(:,find(time>-1.5,1,'first'))/dE); hold on;
plot(energy, n_sat/dE-n(:,find(time>-2.5,1,'first'))/dE); hold on;
ylabel('Population n per level'); xlabel('Energy (eV)')
ylim([-max(n_sat/dE)*0.1,max(n_sat/dE)*1.1]);
grid on; hold off;

figure(4);
grid on; hold off;
ylabel('Stimulated to spontanoues ratio');
xlabel('Energy (eV)');

fluence = sum(p*delt)*energy(center)*1.62e-19 / L^2*1e14 * 1000; %mJ/cm^2
px2=px/px(end)*fluence; %calibrated pump fluence level

figure(5);
plot([0,px2], [0,Pout], 'o-');
ylabel('Output power'); xlabel('Pump fluence (mJ/cm^2)'); grid on;

figure(6);
fwhm2 = fwhm/M*(energy(end)-energy(1));
plot(px2, fwhm2*1e3, 'o-');
ylim([0,100]);
ylabel('Linewidth (meV)'); xlabel('Pump fluence (mJ/cm^2)'); grid on;
